rm(list = ls(all.names = TRUE))
gc()
set.seed(211917)
options(scipen=999)

packages <-c("tidyverse","foreign","stargazer", "lmtest", "multiwayvcov")

new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

lapply(packages, require, character.only = TRUE)
rm(packages, new.packages)

setwd("PUT YOUR WORKING DIRECTORY HERE")

# load data
data <- read.dta("./data.dta")


#---------------#
# Data cleaning # 
#---------------#

data$q6_years_married[data$q6_years_married == -9] <- 0
data$q31_combined_s <- scale(data$q31_combined)
data$q32_combined_s <- scale(data$q32_combined)

## turn into time series
# sorry for the clunkiness!

outcome <- rbind(as.matrix(data$q32_combined), as.matrix(data$q31_combined))
outcome_s <- rbind(as.matrix(data$q32_combined_s), as.matrix(data$q31_combined_s))
treat <- rbind(as.matrix(rep(0,576)), as.matrix(rep(1,576)))
principle <- rbind(as.matrix(rep(0,576)), as.matrix(rep(1,576))) # new by KS

id <- rbind(as.matrix(seq(1:576)), as.matrix(seq(1:576)))
respondent_male <- rbind(as.matrix(data$respondent_male),as.matrix(data$respondent_male))
q1_age <- rbind(as.matrix(data$q1_age),as.matrix(data$q1_age))
respondent_muslim <- rbind(as.matrix(data$respondent_muslim),as.matrix(data$respondent_muslim))
q5_convert <- rbind(as.matrix(data$q5_convert),as.matrix(data$q5_convert))
q6_married <- rbind(as.matrix(data$q6_married),as.matrix(data$q6_married))
q6_years_married <- rbind(as.matrix(data$q6_years_married),as.matrix(data$q6_years_married))
q7_has_job <- rbind(as.matrix(data$q7_has_job),as.matrix(data$q7_has_job))
q8_has_income <- rbind(as.matrix(data$q8_has_income),as.matrix(data$q8_has_income))
q8_income_last_month <- rbind(as.matrix(data$q8_income_last_month),as.matrix(data$q8_income_last_month))
q3_kamba <- rbind(as.matrix(data$q3_kamba),as.matrix(data$q3_kamba))
q3_kikuyu <- rbind(as.matrix(data$q3_kikuyu),as.matrix(data$q3_kikuyu))
q3_luo <- rbind(as.matrix(data$q3_luo),as.matrix(data$q3_luo))
q3_somali <- rbind(as.matrix(data$q3_somali),as.matrix(data$q3_somali))
enumerator_muslim <- rbind(as.matrix(data$enumerator_muslim),as.matrix(data$enumerator_muslim))
enumerator_male <- rbind(as.matrix(data$enumerator_male),as.matrix(data$enumerator_male))
holy_endorser<- rbind(as.matrix(data$holy_endorser),as.matrix(data$holy_endorser)) 
some_message<- rbind(as.matrix(data$any_endorser),as.matrix(data$any_endorser)) # new by KS
secular_endorser <- rbind(as.matrix(data$secular_endorser),as.matrix(data$secular_endorser))
q27_radicalization <- rbind(as.matrix(data$q27_radicalization),as.matrix(data$q27_radicalization))
q30_combined <- rbind(as.matrix(data$q30_combined),as.matrix(data$q30_combined))
q31_combined_s <- rbind(as.matrix(data$q31_combined_s),as.matrix(data$q31_combined_s))
q32_combined_s <- rbind(as.matrix(data$q32_combined_s),as.matrix(data$q32_combined_s))

treatment <- rbind(as.matrix(data$treatment), as.matrix(data$treatment)) 
dataa <- as.data.frame(cbind(outcome, treat, id, respondent_male, q1_age, respondent_muslim, q5_convert, q6_married, q6_years_married, q7_has_job, q8_has_income, q8_income_last_month, q3_kamba, q3_kikuyu, q3_luo, q3_somali, enumerator_muslim, enumerator_male, holy_endorser, secular_endorser, q27_radicalization, treatment))
colnames(dataa) <- c("outcome", "treat", "id", "respondent_male", "q1_age", "respondent_muslim", "q5_convert", "q6_married", "q6_years_married", "q7_has_job", "q8_has_income", "q8_income_last_month", "q3_kamba", "q3_kikuyu", "q3_luo", "q3_somali", "enumerator_muslim", "enumerator_male", "holy_endorser", "secular_endorser", "q27_radicalization", "treatment")

dataa_control <- dataa[which(dataa$secular_endorser!=1 & dataa$holy_endorser!=1),]
dataa_treat <- dataa[which(dataa$secular_endorser==1 | dataa$holy_endorser==1),]

christian_d <- dataa[which(dataa$respondent_muslim==0),]
muslim_d <- dataa[which(dataa$respondent_muslim==1), ]
killing_d <- dataa[which(dataa$treatment==1),]
suicide_d <- dataa[which(dataa$treatment==2),]


#--------#
# Models #
#--------#

##  Regression combined
a1 <- lm(outcome ~ treat + holy_endorser + secular_endorser + respondent_male + q1_age + respondent_muslim + q5_convert + q6_married + q6_years_married + q7_has_job + q8_has_income + q8_income_last_month + q3_kamba + q3_kikuyu + q3_luo + q3_somali + enumerator_muslim + enumerator_male, data=dataa)
summary(a1)
a1$coefficients[2]

# clustering standard errors
a1_clustered <- coeftest(a1, vcov.= function(y) cluster.vcov(y, ~ id, df_correction = T))
a1_clustered


##  Regression No Killing
a2 <- lm(outcome ~ treat + holy_endorser + secular_endorser + respondent_male + q1_age + respondent_muslim + q5_convert + q6_married + q6_years_married + q7_has_job + q8_has_income + q8_income_last_month + q3_kamba + q3_kikuyu + q3_luo + q3_somali + enumerator_muslim + enumerator_male, data=killing_d)
summary(a2)

# clustering standard errors
a2_clustered <- coeftest(a2, vcov.= function(y) cluster.vcov(y, ~ id, df_correction = T))
a2_clustered


##  Regression No Suicide
a3 <- lm(outcome ~ treat + holy_endorser + secular_endorser + respondent_male + q1_age + respondent_muslim + q5_convert + q6_married + q6_years_married + q7_has_job + q8_has_income + q8_income_last_month + q3_kamba + q3_kikuyu + q3_luo + q3_somali + enumerator_muslim + enumerator_male, data=suicide_d)
summary(a3)

# clustering standard errors
a3_clustered <- coeftest(a3, vcov.= function(y) cluster.vcov(y, ~ id, df_correction = T))
a3_clustered


##  Regression Christian
a4 <- lm(outcome ~ treat + holy_endorser + secular_endorser + respondent_male + q1_age + respondent_muslim + q5_convert + q6_married + q6_years_married + q7_has_job + q8_has_income + q8_income_last_month + q3_kamba + q3_kikuyu + q3_luo + q3_somali + enumerator_muslim + enumerator_male, data=christian_d)

# clustering standard errors
a4_clustered <- coeftest(a4, vcov.= function(y) cluster.vcov(y, ~ id, df_correction = T))
a4_clustered


##  Regression Muslim
a5 <- lm(outcome ~ treat + holy_endorser + secular_endorser + respondent_male + q1_age + respondent_muslim + q5_convert + q6_married + q6_years_married + q7_has_job + q8_has_income + q8_income_last_month + q3_kamba + q3_kikuyu + q3_luo + q3_somali + enumerator_muslim + enumerator_male, data=muslim_d)

# clustering standard errors
a5_clustered <- coeftest(a5, vcov.= function(y) cluster.vcov(y, ~ id, df_correction = T))
a5_clustered


stargazer(a1_clustered,a2_clustered,a3_clustered,a4_clustered,a5_clustered)
